so the compiler is not yet used for R's own base and recommended packages. While working on my slides for the upcoming Rcpp workshop preceding R/Finance 2011, I thought of a nice test example to illustrate the compiler. Last summer and fall, Radford Neal had sent a very nice, long and detailed list of possible patches for R performance improvements to the development list. Some patches were integrated and some more discussion ensued. One strand was on the difference in parsing between normal parens and curly braces. In isolation, these seem too large, but (as I recall) this is due some things 'deep down' in the parser. However, some folks really harped on this topic. And it just doesn't die as a post from last week demonstrated once more. Last year, Christian Robert had whittled it down to a set of functions, and I made a somewhat sarcastic post argueing that I'd rather use Rcpp to get 80-fold speed increase than spend my time argueing over ten percent changes in code that could be made faster so easily. So let us now examine what the compiler package can do for us. The starting point is the same as last year: five variants of computing 1/(1+x) for a scalar x inside an explicito Package compiler is now provided as a standard package. See ?compiler::compile for information on how to use the compiler. This package implements a byte code compiler for R: by default the compiler is not used in this release. See the 'R Installation and Administration Manual' for how to compile the base and recommended packages.
for
loop. Real code would never do it this way as vectorisation comes
to the rescue. But for (Christian's) argument's sake, it is useful to highlight differences in the parser. We once again use
the nice rbenchmark package to run, time and summarise alternatives:
This replicates Christian's result. We find that function> ## cf http://dirk.eddelbuettel.com/blog/2010/09/07#straight_curly_or_compiled > f <- function(n, x=1) for (i in 1:n) x=1/(1+x) > g <- function(n, x=1) for (i in 1:n) x=(1/(1+x)) > h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1) > j <- function(n, x=1) for (i in 1:n) x= 1/ 1+x > k <- function(n, x=1) for (i in 1:n) x=1/ 1+x > ## now load some tools > library(rbenchmark) > ## now run the benchmark > N <- 1e6 > benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), + columns=c("test", "replications", + "elapsed", "relative"), + order="relative", replications=10) test replications elapsed relative 5 k(N, 1) 10 9.764 1.00000 1 f(N, 1) 10 9.998 1.02397 4 j(N, 1) 10 11.019 1.12853 2 g(N, 1) 10 11.822 1.21077 3 h(N, 1) 10 14.560 1.49119
k()
is the fastest using curlies, and that explicit exponentiation in
function h()
is the slowest with a relative penalty of 49%, or an absolute difference of almost five seconds between the 9.7
for the winner and 14.6 for the worst variant. On the other hand, function f()
, the normal way of writing things, does pretty
well.
So what happens when we throw the compiler into the mix? Let's first create compiled variants using the new cmpfun()
function and then try again:
Now things have gotten interesting and substantially faster, for very little cost. Usage is straightforward: take your function and compile it, and reap more than a threefold speed gain. Not bad at all. Also of note, the difference between the different expressions essentially vanishes. The explicit exponentiation is still the loser, but there may be an additional explicit function call involved. So we do see the new compiler as a potentially very useful addition. I am sure more folks will jump on this and run more tests to find clearer corner cases. To finish, we have to of course once more go back to Rcpp for some <emph>real</emph> gains, and some comparison between compiled and byte code compiled code.> ## R 2.13.0 brings this toy > library(compiler) > lf <- cmpfun(f) > lg <- cmpfun(g) > lh <- cmpfun(h) > lj <- cmpfun(j) > lk <- cmpfun(k) > # now run the benchmark > N <- 1e6 > benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), + lf(N,1), lg(N,1), lh(N,1), lj(N,1), lk(N,1), + columns=c("test", "replications", + "elapsed", "relative"), + order="relative", replications=10) test replications elapsed relative 9 lj(N, 1) 10 2.971 1.00000 10 lk(N, 1) 10 2.980 1.00303 6 lf(N, 1) 10 2.998 1.00909 7 lg(N, 1) 10 3.007 1.01212 8 lh(N, 1) 10 4.024 1.35443 1 f(N, 1) 10 9.479 3.19051 5 k(N, 1) 10 9.526 3.20633 4 j(N, 1) 10 10.775 3.62673 2 g(N, 1) 10 11.299 3.80310 3 h(N, 1) 10 14.049 4.72871
Rcpp still shoots the lights out by a factor of 80 (or even almost 120 to the worst manual implementation) relative to interpreted code. Relative to the compiled byte code, the speed difference is about 25-fold. Now, these are of course entirely unrealistic code examples that are in no way, shape or form representative of real R work. Effective speed gains will be smaller for both the (pretty exciting new) compiler package and also for our C++ integration package Rcpp. Before I close, two more public service announcements. First, if you use Ubuntu see this post by Michael on r-sig-debian announcing his implementation of a suggestion of mine: we now have R alpha/beta/rc builds via his Launchpad PPA. Last Friday, I had the current R-rc snapshot of R 2.13.0 on my Ubuntu box only about six hours after I (as Debian maintainer for R) uploaded the underlying new R-rc package build to Debian unstable. This will be nice for testing of upcoming releases. Second, as I mentioned, the Rcpp workshop on April 28 preceding R/Finance 2011 on April 29 and 30 still has a few slots available, as has the conference itself.> ## now with Rcpp and C++ > library(inline) > ## and define our version in C++ > src <- 'int n = as<int>(ns); + double x = as<double>(xs); + for (int i=0; i<n; i++) x=1/(1+x); + return wrap(x); ' > l <- cxxfunction(signature(ns="integer", + xs="numeric"), + body=src, plugin="Rcpp") > ## now run the benchmark again > benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), + l(N,1), + lf(N,1), lg(N,1), lh(N,1), lj(N,1), lk(N,1), + columns=c("test", "replications", + "elapsed", "relative"), + order="relative", replications=10) test replications elapsed relative 6 l(N, 1) 10 0.120 1.0000 11 lk(N, 1) 10 2.961 24.6750 7 lf(N, 1) 10 3.128 26.0667 8 lg(N, 1) 10 3.140 26.1667 10 lj(N, 1) 10 3.161 26.3417 9 lh(N, 1) 10 4.212 35.1000 5 k(N, 1) 10 9.500 79.1667 1 f(N, 1) 10 9.621 80.1750 4 j(N, 1) 10 10.868 90.5667 2 g(N, 1) 10 11.409 95.0750 3 h(N, 1) 10 14.077 117.3083
7 comments Liked this article? Click here. My blog is Flattr-enabled.
The text below went out as a post to the
r-packages list a few days ago, but I thought it would make sense to post it
on the blog too. So with a little html markup...
Summary
Version 0.9.0 of the Rcpp package is now on CRAN and its mirrors. This
release marks another step in the development of the package, and a few key
points are highlighted below. More details are in the NEWS and ChangeLog
files included in the package.
Overview
Rcpp is an R package and associated C++ library that facilitates integration
of C++ code in R packages.
The package features a complete set of C++ classes (Rcpp::IntegerVector
,
Rcpp:NumericVector
, Rcpp::Function
, Rcpp::Environment
, ...) that makes it
easier to manipulate R objects of matching types (integer vectors, functions,
environments, etc ...).
Rcpp takes advantage of C++ language features such as the explicit
constructor / destructor lifecycle of objects to manage garbage collection
automatically and transparently. We believe this is a major improvement over
use of PROTECT
/ UNPROTECT
. When an Rcpp object is created, it protects the
underlying SEXP so that the garbage collector does not attempt to reclaim the
memory. This protection is withdrawn when the object goes out of
scope. Moreover, users generally do not need to manage memory directly (via
calls to new / delete or malloc / free) as this is done by the Rcpp classes
or the corresponding STL containers.
A few key points about Rcpp:
which deploys the sugar 'ifelse' function modeled after the corresponding R function. Another simple example isSEXP foo( SEXP xx, SEXP yy) NumericVector x(xx), y(yy) ; return ifelse( x < y, x*x, -(y*y) ) ;
where use the sugar function 'sapply' to sweep a simple C++ function which operates elementwise across the supplied vector. The Rcpp-sugar vignette describes sugar in more detail. Rcpp modules Rcpp modules are inspired by Boost.Python and make exposing C++ functions or classes to R even easier. A first illustration is provided by this simple C++ code snippetdouble square( double x) return x*x ; SEXP foo( SEXP xx ) NumericVector x(xx) ; return sapply( x, square ) ;
which (after compiling and loading) we can access in R asconst char* hello( const std::string& who ) std::string result( "hello " ) ; result += who ; return result.c_str() ; RCPP_MODULE(yada) using namespace Rcpp ; function( "hello", &hello ) ;
In a similar way, C++ classes can be exposed very easily. Rcpp modules are also described in more detail in their own vignette. Reference Classes R release 2.12.0 introduced Reference Classes. These are formal S4 classes with the corresponding dispatch method, but passed by reference and easy to use. Reference Classes can also be exposed to R by using Rcpp modules. Extension packackages The RcppArmadillo package permits use of the advanced C++ library 'Armadillo, a C++ linear algebra library aiming towards a good balance between speed and ease of use, providing integer, floating point and complex matrices and vectors with lapack / blas support via R. Armadillo uses templates for a delayed evaluation approach is employed (during compile time) to combine several operations into one and reduce (or eliminate) the need for temporaries. Armadillo is useful if C++ has been decided as the language of choice, rather than another language like Matlab or Octave, and aims to be as expressive as the former. Via Rcpp and RcppArmadillo, R users now have easy access to this functionality. Examples are provided in the RcppArmadillo package. The RcppGSL package permits easy use of the GNU Scientific Library (GSL), a collection of numerical routines for scientifc computing. It is particularly useful for C and C++ programs as it provides a standard C interface to a wide range of mathematical routines such as special functions, permutations, combinations, fast fourier transforms, eigensystems, random numbers, quadrature, random distributions, quasi-random sequences, Monte Carlo integration, N-tuples, differential equations, simulated annealing, numerical differentiation, interpolation, series acceleration, Chebyshev approximations, root-finding, discrete Hankel transforms physical constants, basis splines and wavelets. There are over 1000 functions in total with an extensive test suite. The RcppGSL package provides an easy-to-use interface between GSL data structures and R using concepts from Rcpp. The RcppGSL package also contains a vignette with more documentation. Legacy 'classic' API Packages still using code interfacing the initial 'classic' Rcpp API are encouraged to migrate to the new API. Should a code transition not be possible, backwards compatibility is provided by the RcppClassic package released alongside Rcpp 0.9.0. By including RcppClassic.h and building against the RcppClassic package and library, vintage code can remain operational using the classic API. The short vignette in the RcppClassic package has more details. Documentation The package contains a total of eight vignettes the first of which provides a short and succinct introduction to the Rcpp package along with several motivating examples. Links Support Questions about Rcpp should be directed to the Rcpp-devel mailing list https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-develyada <- Module( "yada" ) yada$hello( "world" )
Dirk Eddelbuettel, Romain Francois, Doug Bates and John Chambers
December 2010
x
. We then get to the
meat and potatoes and load our
Rcpp package as well as
inline to define the same little test function in C++. Throw in
rbenchmark which I am becoming increasingly fond of for these little timing tests,
et voila, we have ourselves a horserace:
# Xian's code, using <- for assignments and passing x down f <- function(n, x=1) for (i in 1:n) x=1/(1+x) g <- function(n, x=1) for (i in 1:n) x=(1/(1+x)) h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1) j <- function(n, x=1) for (i in 1:n) x= 1/ 1+x k <- function(n, x=1) for (i in 1:n) x=1/ 1+x # now load some tools library(Rcpp) library(inline) # and define our version in C++ l <- cxxfunction(signature(ns="integer", xs="numeric"), 'int n = as<int>(ns); double x=as<double>(xs); for (int i=0; i<n; i++) x=1/(1+x); return wrap(x); ', plugin="Rcpp") # more tools library(rbenchmark) # now run the benchmark N <- 1e6 benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1), columns=c("test", "replications", "elapsed", "relative"), order="relative", replications=10)And how does it do? Well, glad you asked. On my i7, which the other three cores standing around and watching, we get an eighty-fold increase relative to the best interpreted version:
/tmp$ Rscript xian.R Loading required package: methods test replications elapsed relative 6 l(N, 1) 10 0.122 1.000 5 k(N, 1) 10 9.880 80.984 1 f(N, 1) 10 9.978 81.787 4 j(N, 1) 10 11.293 92.566 2 g(N, 1) 10 12.027 98.582 3 h(N, 1) 10 15.372 126.000 /tmp$So do we really want to spend time arguing about the ten and fifteen percent differences? Moore's law gets you those gains in a couple of weeks anyway. I'd much rather have a conversation about how we can get people speed increases that are orders of magnitude, not fractions. Rcpp is one such tool. Let's get more of them.
That took all of five seconds while my computer was also compiling a particularly resource-hungry C++ package.... Just in case we needed another illustration that it is hard to navigate the riches and wonders that is CRAN...R> library(sudoku) R> s <- readSudoku("/tmp/sudoku.txt") R> s [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 8 0 0 0 0 1 2 0 0 [2,] 0 7 5 0 0 0 0 0 0 [3,] 0 0 0 0 5 0 0 6 4 [4,] 0 0 7 0 0 0 0 0 6 [5,] 9 0 0 7 0 0 0 0 0 [6,] 5 2 0 0 0 9 0 4 7 [7,] 2 3 1 0 0 0 0 0 0 [8,] 0 0 6 0 2 0 1 0 9 [9,] 0 0 0 0 0 0 0 0 0 R> system.time(solveSudoku(s)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 8 4 9 6 7 1 2 5 3 [2,] 6 7 5 2 4 3 9 1 8 [3,] 3 1 2 9 5 8 7 6 4 [4,] 1 8 7 4 3 2 5 9 6 [5,] 9 6 4 7 8 5 3 2 1 [6,] 5 2 3 1 6 9 8 4 7 [7,] 2 3 1 8 9 4 6 7 5 [8,] 4 5 6 3 2 7 1 8 9 [9,] 7 9 8 5 1 6 4 3 2 user system elapsed 5.288 0.004 5.951 R>
A few dud universes can really clutter up your basement. - Neal Stephenson, In The Beginning. . . was the Command LineWhat a fun read. It s about technology, sure, but more about culture. Neal takes a good look at operating systems, why we get emotionally involved with them, and why Windows is still so popular. He does this with a grand detour to Disneyland, and a hefty dose of humor. The above quote was from near the end of the book, where he imagines hackers creating big bangs from the command line. He starts out the book from some anecdotes from the early 1970s, when he had his first computer class in high school. His school didn t have a computer, but they did have a teletype (the physical kind that used paper) with a modem link to some university s system. But time on that system was so expensive that they couldn t just dial in and run things interactively. The teletype had a paper tape device. You d type your commands in advance, and it would punch them out on the tape. Then when you dial in, it would replay the tape at high speed . Neal liked this because the stuff punched out of the tape were, actually, bits in both the literal and the mathematical sense. This, of course, led to a scene at the end of the schoolyear where a classmate dumped the bin of bits on the teacher, and Neal witnessed megabytes falling to the floor. Although the book was written in 1999, and needs an update in some ways, it still speaks with a strong voice today and is now also an interesting look at what computing was like 10 years ago. He had an analogy of car dealerships to operating systems. Microsoft had the big shiny dealership selling station wagons. Their image was all wrapped up in people feeling good about their purchase like they got something for their money. And he said that the Linux folks were selling tanks, illustrated with this exchange:
Hacker with bullhorn: Save your money! Accept one of our free tanks! It is invulnerable, and can drive across rocks and swamps at ninety miles an hour while getting a hundred miles to the gallon! Prospective station wagon buyer: I know what you say is true but er I don t know how to maintain a tank! Bullhorn: You don t know how to maintain a station wagon either! Buyer: But this dealership has mechanics on staff. If something goes wrong with my station wagon, I can take a day off work, bring it here, and pay them to work on it while I sit in the waiting room for hours, listening to elevator music. Bullhorn: But if you accept one of our free tanks we will send volunteers to your house to fix it for free while you sleep! Buyer: Stay away from my house, you freak! Bullhorn: But Buyer: Can t you see that everyone is buying station wagons?That doesn t mean that Stephenson is just a Linux apologetic. He points out that the CLI has its place, and has a true love-hate relationship with the text-based config files (remember XF86Config before the days of automatic modelines? Back when you had to get out a calculator and work some things out with pencil and paper, or else risk burning out your monitor?) He points out that some people want to just have the thing work reasonably well. They don t want control in fact, would gladly give it up if offered something reasonably pretty and reasonably functional. He speaks to running Linux at times:
Sometimes when you finish working with a program and shut it down, you find that it has left behind a series of mild warnings and low-grade error messages in the command-line interface window from which you launched it. As if the software were chatting to you about how it was doing the whole time you were working with it.
Even if the application is imploding like a damaged submarine, it can still usually eke out a little S.O.S. message.Or about booting Linux the first time, and noticing all sorts of cryptic messages on the console:
This is slightly alarming the first time you see it, but completely harmless.
I use emacs, which might be thought of as a thermonuclear word processor. . . Microsoft Word, were devoted to features like mail merge, and the ability to embed feature-length motion pictures in corporate memoranda, were, in the case of emacs, focused with maniacal intensity on the deceptively simple-seeming problem of editing text. If you are a professional writer i.e., if someone else is getting paid to worry about how your words are formatted and printed emacs outshines all other editing software in approximately the same way that the noonday sun does the stars. It is not just bigger and brighter; it simply makes everything else vanish. For page layout and printing you can use TeX: a vast corpus of typesetting lore written in C and also available on the Net for free.I love these vivid descriptions: programs secretly chatting with us, TeX being a corpus of typesetting lore rather than a program. Or how about this one: Unix. . . is not so much a product as it is a painstakingly compiled oral history of the hacker subculture. It is our Gilgamesh epic. Yes, my operating system is an oral history project, thankyouverymuch. The book feels like a weird (but well-executed and well-written) cross between Douglas Adams and Cory Doctorow. Which makes is so indescribably awesome that I can t help but ending this review with a few more quotes.
Because Linux is not commercial because it is, in fact, free, as well as rather difficult to obtain, install, and operate it does not have to maintain any pretensions as to its reliability. Consequently, it is much more reliable.
what really sold me on it [Debian] was its phenomenal bug database (http://www.debian.org/Bugs), which is a sort of interactive Doomsday Book of error, fallibility, and redemption. It is simplicity itself. When had a problem with Debian in early January of 1997, I sent in a message describing the problem to submit@bugs.debian.org. My problem was promptly assigned a bug report number (#6518) and a severity level (the available choices being critical, grave, important, normal, fixed, and wishlist) and forwarded to mailing lists where Debian people hang out.That should be our new slogan for bugs.debian.org: Debian s interactive Doomsday Book of error, fallibility, and redemption.
Unix is hard to learn. The process of learning it is one of multiple small epiphanies. Typically you are just on the verge of inventing some necessary tool or utility when you realize that someone else has already invented it, and built it in, and this explains some odd file or directory or command that you have noticed but never really understood before.I ve been THERE countless times.
Note the obsessive use of abbreviations and avoidance of capital letters; this is a system invented by people to whom repetitive stress disorder is what black lung is to miners. Long names get worn down to three-letter nubbins, like stones smoothed by a river.
It is obvious, to everyone outside of the United States, that our arch-buzzwords, multiculturalism and diversity, are false fronts that are being used (in many cases unwittingly) to conceal a global trend to eradicate cultural differences. The basic tenet of multiculturalism (or honoring diversity or whatever you want to call it) is that people need to stop judging each other-to stop asserting (and, eventually, to stop believing ) that this is right and that is wrong, this true and that false, one thing ugly and another thing beautiful, that God exists and has this or that set of qualities.
The stone tablets bearing the Ten Commandments carved in immutable stone the original command-line interface
Apparently this actually works to some degree, for police in many lands are now complaining that local arrestees are insisting on having their Miranda rights read to them, just like perps in American TV cop shows. When it s explained to them that they are in a different country, where those rights do not exist, they become outraged. Starsky and Hutch reruns, dubbed into diverse languages, may turn out, in the long run, to be a greater force for human rights than the Declaration of Independence.
Unix has always lurked provocatively in the background of the operating system wars, like the Russian Army.Available for free online, or as a 160-page book from Amazon.
class BagOf<T> BagOf(Class<T> restrictionClass);is not satisfiable when T is a generic class itself (since there is no ArraySet<Double>.class syntax, for example). The best work-around I know is to drop the T subclassing restriction for restrictionClass:
class BagOf<T> BagOf(Class<?> restrictionClass);The cost is low (obviously no difference at runtime) - you just don't assert that the developer using your class specifies a restriction class derived from the class T used in the generics. That won't prevent certain programming errors such as this anymore
BagOf<Integer> bar = BagOf<Integer>(Double.class)but these shouldn't be too hard to find/fix anyway.Before submitting too clever suggestions, please make sure you've tested them. For example "if (obj instanceof T)" is not valid java code: since generics are implemented by erasure, T cannot be referenced in runtime statements.P.S. It would obviously be nice if the Java syntax would allow Foo<Bar>.class (which at runtime would be the same as Foo.class, and at complie time would have the result type Class<Foo<Bar>>), but currently it does not for all I know.P.P.S. I'm not looking for "Class<? extends T>", that is a different situation. The difficult case is when T is a Generic itself, not a subclass.Update: JM Ibanez pointed me to Neal Gafter's Super Type Tokens, which apparently are the same as TypeLiteral in Google Guice. Thanks!
Next.